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Abstract 

We have conducted a direct numerical simulation (DNS) study of dilute turbulent par- 
ticulate flow in a vertical plane channel, considering thousands of finite-size rigid particles 
with resolved phase interfaces. The particle diameter corresponds to approximately 11 wall 
units and their terminal Reynolds number is set to 136. The fluid flow with bulk Reynolds 
number 2700 is directed upward, which maintains the particles suspended upon average. Two 
density ratios were simulated, differing by a factor of 4.5. The corresponding Stokes numbers 
of the two flow cases were 0(10) in the near-wall region and 0(1) in the outer flow. We 
have observed the formation of large-scale elongated streak-like structures with streamwise 
dimensions of the order of 8 channel half-widths and cross-stream dimensions of the order 
of one half-width. At the same time, we have found no evidence of signiflcant formation of 
particle clusters, which suggests that the large structures are due to an intrinsic instability of 
the flow, triggered by the presence of the particles. It was found that the mean fluid velocity 
proflle tends towards a concave shape, and the turbulence intensity as well as the normal stress 
anisotropy are strongly increased. The effect of varying the Stokes number while maintaining 
the buoyancy, particle size and volume fraction constant was relatively weak. 



1 Introduction 



The flow of a suspension of solid particles through vertically-oriented channels, pipes or ducts 
occurs in a large number of industrial applications, such as fluidized beds, riser iiows and pneumatic 
transport lines. One of the central questions raised by previous studies of these flows concerns the 
nature of the mutual interaction between the turbulent flow field and the motion of the dispersed 
particles. 

Experimental investigations of vertical wall-bounded flow have revealed that the intensity of 
fluid velocity fluctuations can be substantially modifled by the addition of even a relatively small 
volume fraction of heavy particles (e.g. Tsuji, Morikawa, and Shiomi 1984^ 'Rogers and Eaton 
19901 IKulick, Fessler, and Eatonj [ToM) 



Suzuki, Ikenoya, and Kasagi [2000; S ato, Fukuichi, and] 



Hichida[ 2000| Hadinoto, Jones, Yurteri, and Curtis 2005). Depending on the values of the 
various parameters, either enhancement or attenuation of the turbulence intensity is observed. In 
general the turbulence modulation effect can be attributed to a competition between two opposing 
mechanisms: dissipation in the vicinity of particles, and generation of velocity fluctuations due 
to wake shedding (Hetsroni 1989). In addition to these direct consequences, particles can also 



indirectly affect the overall turbulent kinetic energy budget by interacting with specific coherent 
flow structures, thereby changing the very structure of the turbulent flow fleld. Indeed, it has been 
experimentally observed that small particles exhibit so-called 'preferential concentration', e.g. in 



vertical channel flow (Fessler, Kulick, and Eaton 1994) and homogeneous isotropic turbulence 
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(Aliseda, Cartellier, Hainaux, and Lasheras 2002 Wood, Hwang, and Eaton, 2005). Maxey 



(1987) has shown that the accumulation generahy occurs in regions of low vorticity and high 
strain rate. This feature was later confirmed by various DNS studies using the point-particle 



approach ( Squires and Eaton 1990 Wang and Maxey 1993 Pan and Banerjee 1996 Rouson and 



Eaton 2001 Ahmed and Elghobashi 2001 MarchioU and Soldati 2002) 



In 'two-way coupled' point particle simulations the hydrodynamic force acting upon the par- 
ticles is fed back to the fluid, thereby allowing for the description of the influence of dispersed 
particles upon the carrier phase. Using this methodology, some of the mechanisms of turbulence 
modulation through local accumulations of particles has been elucidated in wall-bounded shear 
flow (Pan and Banerjee 1996), homogeneous shear flow (Ahmed and Elghobashi 2000) and ho- 



mogeneous isotropic turbulence ( Squires and Eaton 1990 



Sundaram and Collins 1999 Ferrante 



and Elghobashi 2003). It should be underlined, however, that most of the aforementioned results 
on 'preferential concentration' effects were obtained for particles whose diameter is smaller than 
the smallest flow scales, and for small values of the Reynolds number of the flow around individual 



particles (see also the review by .Eaton and Fessler 1994). 



A different type of mechanism can lead to the formation of particle agglomerations when the 
particle Reynolds number attains values of 0(100). In this regime, where the point-particle ap- 
proach loses its validity, resolved DNS of the sedimentation of finite-size particles in ambient 



fluid (Kajishima and Takiguchi 2002 Kajishima 2004) has revealed that above a critical particle 



Reynolds number of approximately 300, particles tend to form clusters which are elongated in the 
direction of gravity. These clusters have a collective settling velocity exceeding that of an isolated 
particle, thereby inducing flow perturbations at length scales much larger than the particle diam- 
eter. Kajishima (2004) has also shown that angular particle motion plays an important role in the 
regeneration cycle of the agglomerations. Elongated particle accumulations were also reported by 



Nishino and Matsushita ( 2004 ) in experiments on grid-generated turbulence with particle diame- 
ters corresponding to 2-3 Kolmogorov length scales and particle Reynolds numbers in the range 
of 140-210. The formation of these clusters has been attributed to the effect of wake attraction, 
which is the result of a reduced drag force acting on a sphere located inside another sphere's wake. 
The variation of drag and lift experienced by pairs of spheres in different arrangements has been 



thoroughly investigated while keeping the solid objects fixed (Tsuji, Morikawa, and Terashima 



1982 


Tsuji, Morikawa, and Fujiwara 


and Fan 




1994 


). When considering a 



1985 Kim, Elghobashi, and Sirignano[ 1993) Zhu, Liang 



of gravity after initially being vertically separated, the lower drag of the trailing sphere leads to 
an approach of the particles ('drafting'), a subsequent particle encounter ('kissing') and finally a 
lateral separation ('tumbling'), yielding a side-by-side arrangement (cf. Wu and Manasseh 1998). 



Similar particle behavior due to wake effects has been identified in fluidization experiments in- 



volving many spheres (Fortes, Joseph, and Lundgren 1987) 



Despite considerable progress in the study of dispersed two-phase flows, a full understanding 
of the mechanisms driving the interaction between the fluid and solid phases is still lacking at the 
present date, especially in the regime where the smallest length scales of the turbulent flow are 
comparable to the particle diameter. In part this is due to the difficulties which are encountered in 
both laboratory and numerical experiments involving suspension fiow, leading to a relative scarcity 
of detailed data. Concerning fully-resolved direct numerical simulations, the challenge lies in the 
fact that all relevant flow scales (pertaining to the flow around the particles and the background 
turbulence) need to be accurately represented, while at the same time imposing the appropriate 
boundary conditions at the moving phase interfaces. 

High-resolution numerical studies of flnite-size particles interacting with turbulence have often 
focussed on individual spheres. Kim, Elghobashi, and Sirignano (1998) have analyzed a mobile 



spherical particle in synthetic unsteady flow; Bagchi and Balachandar ( 2003 2004 ) have considered 



a flxed particle swept by a homogeneous-isotropic turbulent field; Burton and Eaton (2005) have 
studied the interaction between a fixed particle and homogeneous-isotropic turbulence; |Zeng,] 
Balachandar, Fischer, and Najjar (2007) have determined the effect of wall-bounded turbulence 



upon the flow around a flxed particle. These investigations have provided highly detailed data on 
the time-dependent hydrodynamic forces acting upon a single particle in unsteady or turbulent 
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flow conditions; in the case of Burton and Eaton (2005) the efl'cct of the presence of the particle 
upon the carrier phase was also evaluated. However, true particle mobility while allowing for 
collective effects could not be simulated with these approaches. 

To our knowledge, the only previous direct numerical simulation study of turbulent verti- 
cal channel flow with finite-size mobile particles has been performed by [Kajishima, Takiguchi 



Hamasaki, and Miyake ( 2001[ ). However, these authors considered only a small number (36) of 



relatively large particles (with a diameter corresponding to 32 wall units) and the angular particle 
motion was neglected. It should be mentioned that Pan and Banerjee ( |1997 ) have conducted 
resolved DNS of turbulent particulate flow in a horizontal open channel. They considered up to 
160 stationary and mobile particles, the particle Reynolds number being 0(10) in the latter case. 
Ten Gate, Derksen, Portella, and Van Den Akker ( 2004j ) have performed resolved DNS including a 



dense suspension of finite-size particles (solid volume fraction up to 10~^) in forced homogcneous- 
isotropic turbulence. The particle Reynolds number in their case was 0(10), precluding significant 
wake effects. 

The aim of the present study is first and foremost to generate a body of detailed data describing 
the structure of the carrier phase as well as the motion of the dispersed phase in turbulent wall- 
bounded fiow. For this purpose we have performed interface-resolved DNS of turbulent flow in a 
vertical plane channel configuration involving several thousand spherical particles and integrating 
the equations of motion over O(IOO) bulk flow time units. The parameters are chosen such that the 
mean particle Reynolds number is 136, which means that significant wake effects can be expected. 
The specific questions which we will address concern on the one hand the occurrence of preferential 
concentration of the solid phase, and on the other hand the possible formation of particle-induced 
flow structures. 

The outline of this article is as follows. In the following section the numerical method employed 
in our simulations is explained. We report the various validation checks and grid convergence tests 
which have been performed, specify the conditions of the simulated flow cases and describe the 
initialization procedure. In §[3] we focus on the Eulerian statistics. The procedure for computing 
the statistical quantities and their convergence is described in detail. We discuss the results for 
Eulerian one-point moments and probability density functions of both phases before turning to 
the Lagrangian statistics in § |4] The spatial structure of the dispersed phase is analyzed in § [5] 
and the occurrence of new coherent flow structures is investigated in § [6j The paper closes with a 
short summary and discussion. 



2 Numerical method 



It has been recognized for some time that DNS of particulate flows can be efficiently performed 
on computational meshes which are fixed in space and time, while solving a single set of equations 



in the entire domain, i.e. including the space occupied by the solid bodies ( 


Hofler and Schwarzer 


2000 


Glowinski, Pan, Hesla, Joseph, and Periaux 


2001 


Kajishima and Takiguchi 


2002 


I. In 



this framework, various techniques have been used by different authors in order to impose the 
constraint of rigid body motion upon the fluid at the particle locations. 

Our present simulations have been carried out with the aid of a variant of the immersed 
boundary technique (Peskin 1972 2002) proposed by Uhlmann (2005a I. This method employs a 
direct forcing approach, where a localized volume force term is added to the momentum equations. 
The additional forcing term is explicitly computed at each time step as a function of the desired 
particle positions and velocities, without recurring to a feed-back procedure; thereby, the stability 
characteristics of the underlying Navier-Stokes solver are maintained in the presence of particles, 
allowing for relatively large time steps. The necessary interpolation of variable values from Eulerian 
grid positions to particle-related Lagrangian positions (and the inverse operation of spreading the 
computed force terms back to the Eulerian grid) are performed by means of the regularized delta 
function approach of Peskin (1972 2002). This procedure yields a smooth temporal variation of 



the hydrodynamic forces acting on individual particles while these are in arbitrary motion with 
respect to the fixed grid. 
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Since particles are free to visit any point in the computational domain and in order to ensure 
that the regularized delta function verifies important identities (such as the conservation of the 
total force and torque during interpolation and spreading) , a Cartesian grid with uniform isotropic 
mesh width Ax — Ay — Az is employed. For reasons of efficiency, forcing is only applied to the 
surface of the spheres, leaving the flow field inside the particles to develop freely. IFadlun, Verzicco^ 



Orlandi, and Mohd-Yusof (2000) have d emonstrated that the external fiow is essentially unchanged 
by this procedure; in Uhlmann ( 2005b I it was confirmed that the impact upon the particle motion 
is indeed small. Figure [IJ shows the distribution of forcing points on the surface of a sphere with 
diameter D j Ax = 12.8 and the corresponding Eulcrian grid in the vicinity. 



The immersed boundary technique of Uhlmann ( 2005a ) is implemented in a standard fractional- 
step method for incompressible flow. The temporal discretization is semi-implicit, based on the 
Crank-Nicholson scheme for the viscous terms and a low-storage three-step Runge-Kutta procedure 



for the non-linear part (Verzicco and Orlandi 1996j). The spatial operators are evaluated by central 
finite-differences on a staggered grid. The temporal and spatial accuracy of this scheme are of 
second order. 

The particle motion is determined by the Runge-Kutta-discretized Newton equations for trans- 
lational and rotational rigid-body motion, which are explicitly coupled to the fluid equations. The 
hydrodynamic forces acting upon a particle are readily obtained by summing the additional vol- 
ume forcing term over all discrete forcing points and multiplying the result by a factor of —1. 
Thereby, the exchange of momentum between the two phases cancels out identically and no spu- 
rious contributions are generated. The analogue procedure is applied for the computation of the 
hydrodynamic torque driving the angular particle motion. 

During the course of a simulation, particles can approach each other closely. However, very 
thin inter-particle films cannot be resolved by a typical grid and therefore the correct build-up of 
repulsive pressure is not captured which in turn can lead to possible partial 'overlap' of the particle 
positions in the numerical computation. In practice, we use the artificial repulsion potential of 



Glowinski, Pan, Hesla, and Joseph ( 1999 1, relying upon a short-range repulsion force (with a range 
of 2Aa;), in order to prevent such non-physical situations. Essentially the same method is used 
for the treatment of particles approaching solid walls. The stiffness parameter appearing in the 
above mentioned authors' definition of the repulsion force has been set to ep = 8 • 10^*1?/ (p/Uc oo) 
(where Mc,oo is the terminal particle settling velocity in ambient fluid), based on calibration in 
simulations of two sedimenting particles and particle-wall interactions. It should be noted that 
the use of an approximate treatment of particle 'collisions' somehow taints an otherwise 'direct' 
numerical simulation method. It is also true that the details of individual particle encounters 
sometimes depend on the method for the collision treatment, e.g. in the highly unstable case of 
'drafting-kissing-tumbling' of two particles (Glowinski et al. 2001). However, when considering 



dilute suspensions particle collisions are expected to occur so infrequently that representing their 
exact dynamics does not seem vital for a faithful description of the flow. Previous simulations 



incorporating the method of Glowinski et al. ( 1999 1 have demonstrated that a very good agreement 
with experimental measurements can be achieved even in a relatively dense system ( Pan, Joseph, 



Bai, Glowinski, and Sarin 2002 1. 



2.1 Validation and grid convergence 



The present numerical method has been submitted to exhaustive validation tests ( Uhlmann 2004 
2005a|b 2006a). These tests have established the following features: (i) second-order accuracy 
of the interpolation scheme and low sensitivity of the error with respect to the position of the 
immersed object relative to the grid; (ii) favorable smoothness properties of the scheme in com- 
putations involving oscillating objects. 

In addition, we have performed a detailed grid convergence study of the motion of particles 
in pressure-driven vertical plane channel flow ( Uhlmannl 2006b). Two regimes were considered: 



laminar flow with a single heavy particle of diameter D/h = 1/20 (h being the channel half 
width) at i?e £3^00=100 and turbulent fiow with four identical particles at i?e£)_oo=136 (the terminal 
Reynolds number fle^ oo is based upon the particle diameter and the terminal settling velocity in 



4 



ambient fluid). In the former case, four different grids were considered, corresponding to a particle 
resolution of D//S.X — {12.8, 19.2, 25.6,38.4}, while the CFL number was kept constant. Second 
order convergence of the particle velocities was observed in this test. In the turbulent case, the flow 
conditions were chosen equal to those of the configuration of interest in the following section (cf. 
case B of table [T] and [2]). Three spatial resolutions, D/ /S.x = {12.8, 19.2, 25.6}, were analyzed in a 
computational domain which was sufficiently large to allow for sustenance of the turbulent state. 
During the observation interval of approximately 3.6 bulk fiow time units, all three grids yield the 
same flow features. The instantaneous particle velocities obtained from the coarsest and the finest 
grid match to within 8.3% of the terminal velocity; the second-moment fluid statistics exhibit 
a maximum relative discrepancy of 7%. It should be noted that particle paths computed with 
two different spatial resolutions will invariably deviate with time (even when both simulations are 
extremely well resolved) since perturbations can be amplifled under supercritical (i.e. turbulent) 
conditions. From these considerations and the above mentioned comparison with experimental 
data for sedimentation of a single sphere, it was concluded that a resolution of D/ Ax = 12.8 is 
sufficient for our present investigation of turbulent particulate channel flow. Figure [T] shows this 
grid in the vicinity of a sphere. 

In separate simulations it was verifled that our DNS code reproduces single-phase turbulence 
results for plane channel flow faithfully. In particular, using a mesh width of Ax/h — 1/256 (which 
corresponds to the value chosen in the following section) gives results in very good agreement with 
the reference data of Kim, Moin, and Moserj (1987). The mean velocity profile from our single- 
phase computation has been included in figurel6F 



2.2 Flow configuration 

We are considering particulate flow in a plane channel which is aligned along the gravitational (x) 
direction {y is the wall-normal and z the spanwise coordinate, cf. figure [2]). The fiuid is driven 
upwards by a mean pressure gradient {p).x < 0. The computational domain is periodic in the 
streamwise and spanwise directions (with periods and L^), i.e. particles leaving the domain on 
any of the four periodic end-planes re-enter at the opposite side. At the solid walls, the no-slip 
condition is imposed upon the fluid, while the particles rebound due to the short-range repulsion 
force, as specified in S|2j The bulk flow Reynolds number is maintained constant at a value of 
Rcb = Uf, h/v = 2700 (uf, being the bulk velocity and v the kinematic viscosity), which generates 
a turbulent flow with Rcr — Ur h/v 172 [ur being the wall shear velocity) in the absence of 
particles. The particle diameter is chosen as D/h ^ 1/20, corresponding to D'^ = Dur/v = 8.6 
in wall units of the unladen flow (11.25 wall units of the laden flow). Physically, this regime of 
particle sizes is interesting, since the particle diameter is comparable to the cross-stream length 
scales of coherent structures typically found in the buffer layer of wall-bounded turbulent flow; a 
mutual interaction can therefore be expected. Numerically, this particle size is the smallest value 
for which we can accurately represent the local flow at a particle Reynolds number of 0(100) on 
the grid required for an accurate computation of the single-phase turbulence. As mentioned in 
the previous section, we choose a spatial resolution of D/ Ax — 12.8, which corresponds to a mesh 
width of Ax/h — 1/256, i.e. Ax'^ = 0.67 in wall units. The time step for case A (the two flow 
cases are further specified below) was set to At = 0.5761Aa;/Mft, which corresponds to a maximum 
CFL number of approximately 1; for case B the step was reduced to At = 0.4148Aa;/ub, i.e. a 
maximum CFL number of approximately 0.75. 

We have set the bulk fluid velocity {ub) equal to the terminal particle velocity (uc,oo)- This 
condition means that the average vertical velocity of the particles will be close to zero, supposing 
that their average wall-normal spatial distribution is near uniform and neglecting wall effects, 
collective effects and possible effects arising from the interaction with turbulence. With these 
settings it turns out that the vertical velocity averaged over all particles takes a slightly positive 
value as we will see below. From the condition Uc.oc^ui, it follows that the terminal particle 
Reynolds number measures Reo.oo ~ 136. From the equality between buoyancy and drag forces 
acting upon an isolated particle, we can form the Archimedes number, which takes the following 
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= (D//i)^i?eg^ ( ^ - 1 )= 13328, (1) 



value: 

Ar = (D/hf Rel ^ ( 

Pf 

where Pp/ Pf is the density ratio between particle and fluid and g — {gx, 0, 0) the vector of gravita- 
tional acceleration (with Qx < 0). This leaves us with one parameter free to choose {pp/ Pf^ say), 
the other one {\gx\h/u1) being fixed by ([l]). 

We have analyzed two cases with parameter values as given in tabic [l] The current range 
of density ratios is representative of the motion of glass, ceramics or sand in water and is of- 
ten investigated in laboratory experiments (e.g. Parthasarathy and Faeth[ |1990 Suzuki et al 



2000 Kiger and Pan 2002). Apart from its relevance for the study of fundamental questions 



in turbulence-particle interaction, this regime has important applications in environmental flows 
and hydrology. The relative particle density in case A is approximately 4.5 times higher than 
in case B, and the Stokes number (defined as the ratio between particle and fluid time scales) 
varies accordingly. In particular, we determine the particle time scale from the usual Stokes drag 
law, TsD = -D^pp/(18z^p/), and consider two time scales for the fluid motion: the near-wall scale 
T+ — v/u^ and the bulk flow scale Tf, = h/ub- Consequently, we can form two Stokes numbers, 
characterizing the ratio of time scales in the near- wall region, St~^ = tsd/t"*", and in the outer 
flow Stb = Tso/Tb- Table [T] shows that the present particles have a relatively large response time 
with respect to the flow in the near-wall region, where 5*^+ > 10, but the scales of both phases 
are comparable in the bulk of the flow, i.e. Stb ~ C(l). 

The global solid volume fraction has been set to = 0.0042, which can be considered at the 
upper limit of the dilute regime. 

The two cases were run in computational domains of different wall-parallel dimensions, the 
volume in case B being eight times larger {Ah x 2h x Ih versus 8h x 2h x 4/i). The number of 
particles in cases A, B is therefore Np — {512,4096}, respectively. The domain used in case A 
evidently constitutes a smaller sample size per instantaneous flow field. However, it allows for 
simulations over longer time intervals which is of special interest for the computation of slowly 
evolving statistical quantities acquired along the particle paths. Table [2] shows the global grid size 
of the simulations, which exceeds 10® Eulerian nodes and 2 • 10^ Lagrangian force points in case B, 
where the work is typically distributed over 512 PowerPC processors with fiberoptical interconnect 
(MareNostrum at BSC-CNS, Barcelona). It should be mentioned that the simulation of case B 
has required roughly 10^ CPU hours of execution time. 



2.3 Initialization of the simulations 

Particles are introduced into fully-developed turbulent single-phase flow at time Iq. The initial 
particle positions form a regular array covering the computational domain, perturbed randomly 
with an amplitude of one particle diameter. Their initial velocities are matched with the fluid 
velocities found at to at the locations of each particle's center; the initial angular particle velocities 
are set to zero. During the first time step of the two-phase computation, the coupling algorithm 
solidifies the fluid occupying the volume of each particle and buoyancy forces set in. Figure [sj^a) 
shows how the average over all particles of the vertical particle velocity drops from a value close 
to the bulk fluid velocity to almost zero over an initial interval of approximately 3 (1) bulk time 
units in case A (B). Let us mention that the long-time averaged mean upward drift of the particles 
amounts to 0.029u(, and 0.058Mb, respectively. The data corresponding to the initial transient will 
not be considered in the following and, therefore, an initial interval of 17 bulk time units has been 
discarded. The length of the observation interval for the two cases is indicated in table [2j 

From figure[3j[6) it can be seen that the friction velocity increases immediately after the addition 
of particles, leading to a sharp rise of i?e^ which then oscillates around a new mean value of 
approximately 225 in both cases. The amplitude and frequency of the oscillations is larger in case 
A due to the smaller box size. This higher friction velocity in the particulate cases rcfiects the 
modification of the mean fluid velocity profile as discussed below (cf . figure [6]) . From here on, 
it is understood that quantities given in wall units (indicated by the usual superscript '+') are 
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normalized with the actual friction velocity, i.e. the increased values are used in the particulate 
flow cases. 



3 Eulerian statistics 

3.1 Computation and convergence of statistical quantities 

Let us introduce the notation for the different averaging operations employed in the following. 
(Ox refers to the streamwise average of a quantity ^, {^)xz is the average over wall-parallel planes 
and (^)t is the temporal average (over the observation interval given in table[2]); the average over 
time and wall-parallel planes is denoted by (^) — {{(,)xz)t', finally, (^)p refers to the instantaneous 
ensemble average over all particles. 

The Eulerian space averages of quantities related to the dispersed phase are carried out for 
discrete bins in the wall-normal direction, defined as lyj = [j — ■ 2h/Ni,in VI < j < -/Vbhi- In 
practice a value of N^in = 80 was chosen (except where otherwise stated), such that the width of 
each bin corresponds to one particle radius. 

The time- and plane-averaged number density is then computed for each bin by the following 
formula: 

, nabs 

(-^)(y^) = ^ , L L \i ^ ^' 

where Xc*'' = (a;!*"* , j/c*'' , a^c*"* ) are the center coordinates of the zth particle, yj is the coordinate of 
the center of the bin lyj , Uobs the number of available records during the observation interval and 
\Iyj \ is the width of jth wall- normal interval. The solid volume fraction simply becomes: 

(0,)(y,) = K(n,)(y,), (3) 

denoting the particle volume by Vc = D^tt/6. 

The time- and plane- wise average of a quantity is computed as follows: 

fe>"'''- (,..>.„.i„...z..L.i^.i i: ^ 



1=1 i=l 

yi'''{t,)<£lyj 



Here can stand for any component of the linear or angular particle velocity vectors, which are 
denoted by ui*^ and u;c \ respectively. 

In order to analyze the motion of the dispersed phase, we have stored the full set of particle- 
related quantities after every 10 (5) time steps for case A (B); this corresponds to approximately 
1.5 • 10'' (5.2 • 10^) samples. Statistical quantities were then computed in a post-processing stage. 

In the vicinity of sharp gradients, discrete binning of the data for the dispersed phase tends 
to smooth out the profiles. This effect can be seen in figure Qa), where the solid volume fraction 
of case B is shown for different numbers of wall-normal bins varying from iVf,i„ = 40 to 320. The 
plot confirms the adequacy of the overall number of discrete samples, because the profile of {4>s) 
does not become significantly more noisy when increasing the number of bins. Since the particles 
have a non- negligible size with respect to the channel width [D/h — 1/20), the apparent drop 
from the maximum to a value of zero at the wall could be an artifact of the binning. Figure |4j6) 
demonstrates that this is not the case. Using fine bins with a width of D/8 reveals first that there 
are no particles within yc < -D/2, which is the zone 'forbidden' by the geometric constraint of the 
wall. Secondly it shows that the solid volume fraction smoothly increases for wall-distances larger 
than this lower limit. Furthermore, let us point out that the mean wall-normal particle velocity is 
neghgible (max | (i;c) |/'"b < {10~^, 7 • 10~^} in cases A, B), indicating that the particle distribution 
has indeed achieved a statistically steady state. 
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Concerning the statistics of the carrier phase, we have accumulated the usual low-order one- 
point moments (mean velocity, stress tensor) during the simulations. For this purpose, averaging 
was performed over the entire composite flow field, i.e. containing the regions of the immersed solid 
particles (the average velocity obtained by this procedure is denoted by (u) etc) . This decision was 
taken in favor of computational efficiency, and its consequences are discussed below. On the other 
hand, we have stored a number of full instantaneous flow flelds and corresponding particle data 
for visualization purposes and a posteriori computation of additional quantities, such as two-point 
correlation functions. From these flelds wc have estimated the difference between averaging flrst 
and second moments over the composite flcld versus taking into account the actual fluid nodes 
only; the fluid-only averaged velocity is denoted by {u)f. For more details the reader is referred 
to appendix [A] This test has lead us to the conclusion that the employed method overestimates 
the streamwise velocity fluctuations (by up to 6.7%) and slightly underestimates the mean fluid 
velocity (by less than 0.7%); the wall-normal fluctuation and the Reynolds stress are overestimated 
by approximately 2%, the spanwise component by 0.7%. These differences affect the data presented 



in some of the figures discussed in the present section as well as the following section § 3.2 (i.e. the 
figures [5j [6) [9) [10|), [iT] [l2j [13}}). As shown in appendix |A) the over-prediction of the streamwise 



velocity fluctuation amplitude can be quite accurately represented by a linear function of the solid 
volume fraction, allowing for a correction of the run-time statistics. This correction has been 
included in flgures 11 12 and [13} b). 

The usual check of the statistical convergence of low-order moments in single-phase plane 
channel flow includes a veriflcation of the linearity of the total shear stress profile. In the multi- 
phase case, however, the streamwise momentum balance includes an additional contribution. The 
equation, as derived in appendix |B| reads: 

dy Rcr dy^ \pf J 

where y = y/h, and (u',w',w') are the components of the fluctuations with respect to the mean 
velocity ((m),0,0). It should be noted that this equation is written for the composite flow fleld 
(including the fluid located at the positions of the solid particles), and all averaging operators 
extend over the entire computational domain fl. It can be seen from ([sj) that whenever the 
particle distribution is not uniform (i.e. when $s — is different from zero), the total shear 
stress does not vary linearly with the wall-distance. In order to test whether the accumulated 
data has indeed achieved statistical equilibrium, we have compared the terms of the integral of 



equation (|5|, viz. 



1 d{u)+ ( Pp \ 9xh 



T tot Tj 

In practice, the integral of the solid volume fraction was evaluated numerically (after piecewise 
polynomial interpolation of the statistical data, using cubic spline functions). Figure [5|^a-6) shows 
the different contributions to the balance (|6|. It can be seen that the profile of Ttot deviates 
quite significantly from the usual straight line, the particle contribution being significant at all 
wall-distances. However, the sum of the stress terms and the particle contribution T^ot + Irf> does 
indeed follow the line y — 1 quite closely. It is further observed that the data in figure [5ja-fo) does 
not verify the odd symmetry with respect to the center at y = 1 as much as would be expected 
from the length of the respective observation intervals and the size of the computational boxes. 
We will return to this lack of symmetry in § J6| below. Additional averaging over both channel 
halves clearly improves the match (cf. figure [5p-(i) . Due to the larger sample size (both for the 
carrier phase as well as the dispersed phase), the statistical data for case B is found to be closer 
to statistical equilibrium than case A. Therefore, in the following we will mainly discuss results 
from case B. It should be noted that in both cases, a certain deviation from linearity is observed 
in a limited region near the wall {y < 0.1). This discrepancy can be attributed to the effect of 
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discrete binning of the solid volume fraction, as discussed above (cf. figure |4p). Finally, it can be 
concluded that the statistical data — when averaged over both channel halves — has approximately 
achieved a steady state consistent with equation m«. 



3.2 One-point moments 

The mean fluid velocity profiles (computed by averaging over the composite flow field) normalized 
by the bulk velocity are shown in figure [6f^a). We observe a clear difference with respect to the 
single-phase case, characterized by a tendency of the particulate flow to form a concave profile with 
higher gradients at the walls and a flat section near the center of the channel. The concave shape 
is slightly more pronounced in the high-inertia case A. On the other hand, when normalized in wall 
units (cf. figure |6f)), the particulate flow cases exhibit a considerably lower mean velocity than the 
single-phase counterpart across the whole channel. As expressed by the streamwise momentum 
balance equation ([s]) the mean velocity of the composite field is determined by the Reynolds shear 
stress and the particle distribution. Before turning to the discussion of these two quantities, let 
us mention that concave mean velocity profiles have previously been observed in experiments on 
dilute particulate flow at various parameter points: in vertical pipes (|Tsuji et al. 1984 Hadinoto 



et al. 



2005) and — to a lesser extent — in vertical channels (Sato et al. 2000). However, neither the 
distribution of the solid volume fraction nor the Reynolds shear stress profiles are available from 
these experimental data-sets. 

The profiles of the mean solid volume fraction are shown in figure [7j In both cases we observe 
a sharp peak at a wall distance of approximately 16 wall units {y = 0.07), a slightly below-average 
value around y = 0.2 and finally a mild local maximum at the center of the channel. A sharp rise of 



the solid volume fraction near the wall was also reported by Suzuki, Ikenoya, and Kasagi ( 2000 ) in 
their experiment on downward channel flow containing particles with a diameter of approximately 
4 wall units. 

A number of mechanisms which influence the wall-normal distribution of particles have been 
discussed in the literature. One of the principal contributions has been termed turbophoresis 
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opposite to gradients in the turbulence intensity. For an individual particle this means that it is less 
probable to receive the necessary momentum driving it from a region of low turbulence intensity 
towards a high-intensity region than vice versa; turbophoresis is then a statistical consequence. An 
expression corresponding to a turbophoretic force appears formally when considering a Eulerian 
formulation for the dispersed phase. In this context the wall-normal particle momentum equation 
contains the following term: —d{v'^v'^ / dy ( [Young and Leeming 19971. The profile of this quantity 
is shown in figure [s] (the correlation {v'^v'^ itself is shown in figure 13 and will be discussed below). 
It can be seen that the turbophoresis term is positive for wall-distances below « 16 (i.e. driving 
particles towards the center) and then negative (directed wall-ward) up to approximately y = 0.2; 
for larger wall- distances its values are close to zero. The point where the turbophoretic force 
changes sign corresponds to the near-wall maximum of the mean solid volume fraction (^s) (cf. 
figure [t]), and the interval over which the force has significant values matches the region where the 
corresponding bump in (</)s) is observed. This close correspondence suggests that turbophoresis is 
indeed important for the mean wall-normal particle distribution in the presently investigated cases. 
Other effects — such as the mean shear-induced lift force and the turbulent diffusion — are expected 
to contribute as well. However, the additional mechanisms are not directly quantifiable from the 
available data, and, therefore, it is not possible to further investigate their relative importance in 
establishing the equilibrium distribution. 

The Reynolds stress normalized by the friction velocity is shown in figure ^a) . The proflles of 
both simulated cases are of a much smaller amplitude than the single-phase reference data across 
most of the channel width. When normalized by the bulk velocity, the reduction is limited to wall- 
distances y ^ 0.2. Notably, near the center of the channel a considerable region with vanishing 
Reynolds stress is observed. This means that the correlation between the streamwise and the wall- 
normal fluid velocity fluctuations is signiflcantly reduced in the present cases. However, the sum 
of the Reynolds shear stress and the contribution from the solid volume fraction to the integrated 
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momentum balance (|6]), i.e. {u'v')^ + I^, (shown in figure [of)) is of larger absolute value than the 
Reynolds shear stress of the single-phase flow when normalized in wall units. From equation (|6]) it 
follows then that the wall- normal gradient of the mean streamwise velocity (in wall units) is smaller 
in the particulate case, leading to the reduced value of {u)^ observed in figure [6]j6). Furthermore, 
we note from figure [9](6) that the largest difference between the sum -f and the Reynolds 

shear stress of the single-phase case is obtained around y = 0.075, where the profile of the former 
quantity exhibits a visible bump. This is evidently the signature of the near-wall peak of the solid 
volume fraction (cf. figure [?]) which is found at the same location. Consequently, the accumulation 
of particles in the buffer region has the effect of bringing more momentum towards the wall. This 
in turn causes the flattening of the mean velocity profile and generates a shear-free region near 
the center of the channel. 

Figure 10 a) shows the profiles of the mean (streamwise) particle velocity, which is positive in 
the center of the channel and changes sign at approximately y = 0.2. In both cases the shape of 
the profile is very similar to the corresponding mean fluid velocity proflle (cf. figure [sji), but with 
a negative shift. Close to the wall, the mean particle velocity profile differs in shape from the 
fluid velocity counterpart, as will be discussed below. In figure 10 '6) the difference between the 
mean velocities of the two phases, (uc) — (u), is plotted; it is found to be approximately constant 
and nearly equal to —U},. For wall-distances y < 0.1, however, the magnitude of the relative 
mean velocity substantially decreases. This result indicates that the coefficient of the drag force 
acting on average upon the particles in this region is higher than at large wall-distances. Similar 
behavior is reported by Zeng, Balachandar, and Fischer (2005) who have investigated the drag 
force on a single sphere translating in a channel with quiescent fluid. These authors have detected 
a significant increase of the drag coefficient for wall distances below y/D = 2 (i.e. for y < 0.1 
in the present context), at similar values for the particle Reynolds number. A drag coefficient 
based upon the relative mean velocity of our simulations can be computed by formulating the 
equilibrium between buoyancy and drag forces, Cn.present = {Pp/Pf — ^)\gx\^D/{{uc) — (w))^, and 
substituting the data of figure 10 [b). The ratio between this relative mean DNS drag coefficient and 



the standard drag coefficient of an isolated sphere given by the Schiller-Naumann formula (Clift, 



Grace, and Weber 1978) is Cd, present /Cd,sn = 1.46 at a wall distance of y/Z) = 0.9375 (where the 



particle Reynolds number based upon the relative mean velocity is 107.2). As a comparison, the 
value obtained from the correlation formula given by Zeng et al. (2005) is Cd,zbf/Cd,sn — 1-44. 
The match is indeed remarkable, indicating that results obtained for the wall-effect on individual 
particles at ambient ffow conditions are indeed relevant to turbulent particulate flow in the dilute 
regime. 

The r.m.s. fluctuations of fluid velocities, normalized by the friction velocity, are shown in 
figure |11[ This figure also includes profiles of the streamwise component from which the overes- 
timation due to the computation of statistics over the composite flow field has been subtracted 
(using the fit of appendix |A]). The streamwise component in both particulate flow cases is strongly 
enhanced as compared to the single-phase flow, except for a small interval very close to the wall. 
The wall- normal component is slightly increased for y < 0.1 and slightly decreased for wall- 
distances above that value; the spanwise component is decreased almost across the whole channel 
width. When normalized by the bulk velocity, all three proflles exceed the single-phase reference 
data, except for the spanwise component for intermediate wall-distances (0.1 < y < 0.6). Since 
the streamwise component is dominant, the turbulent kinetic energy is significantly enhanced with 
respect to the single-phase flow, as shown in figure [T2j When subtracting the overestimation due 
to the composite flow statistics, the ratio reaches approximately 4.6 (5.3) on the centerline in case 
A (B); when normalized in bulk units, the ratio is 7.9 (9.0). The principal reason for the observed 
turbulence enhancement is the generation of velocity fluctuations in the wakes trailing the parti- 
cles. Wakes are indeed found to be the most prominent flow structures, as can be concluded from 
flow visualisations (cf . figure [26] below) . The wake structures significantly increase the turbulence 
intensity, but they do not contribute to the generation of the mean Reynolds shear stress, since 
correlations between u' and v' internal to each wake cancel out upon averaging over an ensemble 
of mostly uncorrelated wakes (the random-like distribution of particles will indeed be confirmed in 
§[5|. In addition to the wakes, new large-scale velocity structures representing streamwise velocity 
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perturbations appear in the particulate flow. These flow structures, which will be discussed in §[6j 
further contribute to the enhancement of turbulence intensity. In case A somewhat lower values for 
the streamwise turbulence intensity are obtained when compared to case B. This difference might 
be a consequence of the smaller domain size in the former case, where the large particle-induced 
structures are strongly constrained. A significant increase in streamwise velocity fluctuation levels 
has also been detected by Suzuki et al. ( 2000) , although their solid volume fraction was one order 
of magnitude lower. The authors already attributed this observation to the presence of the wakes. 
The r.m.s. velocity fluctuations of the dispersed phase are shown in figure 13 |a); for conve- 



nience, a comparison of the r.m.s. velocity fluctuations of both phases in case B (i.e. combining 
data of figures 11 and 13 1) is shown in figure 13 6). A similarity between the corresponding com- 
ponents of the two phases can be observed on these graphs, except for the immediate vicinity of 
the wall. The streamwise particle r.m.s. velocities also exhibit a mild peak in the buffer layer, and 
the shape of the wall-normal and spanwise components is similar as for the fluid data. This simi- 
larity shows that the present particles are able to respond to the fluid velocity fluctuations. The 
kinetic energy of the particle velocity fluctuations (not shown) is lower than the fluid turbulent 
kinetic energy across the whole channel. This behavior is expected due to the particle inertia. 
In the high-inertia case A, the intensity of the fluctuations of all particle velocity components is 
smaller than the fluid data. In case B, the streamwise component is also considerably smaller 
than the corresponding fluid values. The level of the remaining two components, however, slightly 
exceeds the fluid counterparts across most of the channel. In addition, it can be seen in flgure [Tsl 
that the particle velocity fluctuation intensities flatten when approaching the wall and increase 
again for small wall-distances y ^ 0.06 (corresponding to 14 wall units). What looks like a 'kink' 
of the curves is not an artifact of the data analysis, since a reflnement of the width of the bins 
used in the computation of the Eulerian statistics has conflrmed the observed trend. It should be 
remarked that the particles are not subject to the no-slip condition at the wall, and, therefore, the 
asymptotic behavior of the particle velocity is expected to differ from the one of the fluid velocity 
(cf. Kulick et al. 1994). However, it is difficult to imagine how the intensity of the fluctuation 
components can rise near the wall, where the flow is increasingly quiescent. Also, the observed 
behavior does not seem to be due to the numerical treatment of particle- wall collisions, since those 
only concern the wall-normal component of the particle motion, whereas the profiles of all com- 
ponents exhibit a similar upward trend. It must be concluded that the reasons for the increased 
near-wall variance of particle motion remain unclear at the time being, requiring further research. 

The profiles of the spanwise component of the mean angular particle velocity (wc,z) are shown 
in figure 14 a). In both cases there is a minimum at a wall-distance of approximately 25 wall 
units, and the profile tends towards zero at the centerline and at the wall. In figure 14^6) the 
DNS data are compared to the mean shear profile. In a linear shear flow the steady rotation of a 
cylinder is described by the formula —Ad{u) /dy, where A is positive and decreases with the shear 
Reynolds number (cf. Ding and Aidun 20001 Zettner and Yoda 2001 ). Constrained simulations of 



a single neutrally-buoyant and spherical particle placed in laminar pipe flow have likewise yielded 
a proportionality between the steady-state angular particle velocity and the wall-normal fluid 
velocity gradient over most of the pipe radius (Yang et al. 2005[ figure 3). It can be seen In 
figure 14 5) that our data for case B are reasonably well represented when choosing A = 0.1, 
except for the region y'^ < 25 (similarly for case A, which is not shown). The discrepancy in the 
near-wall region can probably be attributed to a modification of the fiow field around particles 
in the immediate vicinity of the wall induced by the no-slip condition. A similar tendency was 



observed for a cylinder in linear shear fiow when the confinement ratio was lowered (Ding and 



Aidun 2000) 



3.3 Probability density functions and quadrant analysis 

Probability density functions (pdfs) of the fluid velocity components were computed from 12 
instantaneous flow flelds of case B, taking into account only values at grid nodes inside the fluid 
domain. The sample size was doubled by making use of the symmetry with respect to the center of 
the channel. In the following we will discuss the normalized pdfs in the buffer region (y+ — 20) and 
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at the center of the channel, in comparison to the single-phase data of |Moser, Kim, and Mansour| 
( 1999 ) at Rct = 180. The values of the skewness at these wall-distances are given in table [3j 
It can be seen that the skewness of the spanwise velocity component computed from our data 
set is close to zero as a consequence of the spanwise periodicity, indicating that sufficient sample 
has been obtained. In figures 15 and 16 we observe that the pdf for the streaniwise component is 
strongly non-symmetric with respect to the mean value. While the branch for positive fluctuations 
matches those of the single-phase data closely, there exists a higher probability of finding strongly 
negative fluctuations in the present particulate case. This tendency is confirmed for both wall- 
distances and it reflects the negative perturbations induced by the streamwise momentum deficit 
in the wakes trailing the particles. The comparison with the single-phase data shows that the 
infiuence of the particles upon the shape of the pdfs for the wall-normal and spanwise velocity 
components is relatively small, manifesting itself through a small increase of the non-Gaussian 
tails. The skewness of the wall-normal component in the buffer layer is reduced from its value in 
single-phase flow. 

The corresponding pdfs of the particle velocity are shown in figure[T7]for different wall-distances 
(in the buffer region, the logarithmic layer and in the center of the channel) ; the skewness values are 
given m table m The curves for the streamwise component are found to largely resemble Gaussian 
functions. However, they exhibit discernible non-Gaussian tails on the negative side which grow 
with the wall-distance. This feature indicates that the particles respond to some degree to the 
negative fluid velocity fluctuations which they experience in the wakes of other particles. On the 
other hand, the pdfs for the wall-normal component are characterized by much more pronounced 
non-Gaussian tails. Additionally, in the buffer region one can observe a strong skewness of v'^, 
towards negative values, i.e. high-velocity motion (with intensity above four standard deviations) 
directed towards the wall is occurring more frequently than away from it. Since the pdf of wall- 
normal fluid velocity at the same wall-distance was not found to exhibit a comparable skewness 
(cf. figure 15), comparison between tables [s] and |4]) , the present result suggests that the particles 
respond selectively to certain coherent structures in the buffer layer. Data accumulated for the 
wall-normal particle velocities in case A (not shown) features a less pronounced negative skewness. 
Unfortunately, no experimental data for the skewness of the wall-normal particle velocity pdf in 
channel flow could be obtained. Therefore, in order to clarify this point, more in-depth investiga- 
tions of the interaction between buffer-layer structures and particles with different values of the 
Stokes number should be undertaken. In particular, conditional averaging of wall-normal particle 
motion based upon the type of the surrounding flow structures (similar to the point-particle study 
of Marchioli and Soldati 2002) should provide valuable insight. Finally, we observe in figure 17 'c) 
that the pdf of the spanwise component of the particle velocity (having a skewness near zero by 
virtue of symmetry) takes a shape very similar to the one found for the wall-normal component 
outside the buffer layer. 

In order to analyze the different contributions to the Reynolds shear stress, we have computed 
joint probability density functions of the fluid velocity fluctuations and v'j:. The data for 
wall-distances of y+ — 20 and 74 is shown in figures 18 a, 6); an exhaustive discussion of the 
corresponding single-phase case can be found in Kim et al. ( 1987[ ). Quadrants are defined in the 



usual manner, such that e.g. the second quadrant corresponds to an ejection of low-speed fiuid (a 
positive value of coupled with negative u^). The joint pdfs for the present case B are found 
to take an elliptical shape, elongated in a direction which is slightly inclined with respect to the 
axis v'j: — 0. Therefore, the contributions from the second and the fourth quadrants dominate 
as in single-phase fiow. It can be seen in both figures |18[ a, b) that the most extreme negative 
fluctuations of are not significantly correlated with fluctuations of any sign, therefore not 
contributing to the mean Reynolds shear stress. In the logarithmic region the joint pdf tends 
towards a less anisotropic shape. 

The fractional contributions from each quadrant to the Reynolds shear stress are given in 
figure 19 In the buffer layer (y+ = 20) and in the logarithmic region (y+ = 74) the total 



contribution from the second quadrant is the strongest, with the fourth quadrant closely following. 
The main difference with respect to the single phase data (as discussed in Kim et aTj 19871 is 
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the relative importance of the quadrants one and three (i.e. those which correspond to positive 
Reynolds shear stress). In the particulate case, the fractional contributions of all four quadrants are 
significantly higher than in the case without particles, most probably due to fluctuations generated 
in the wake of particles. As discussed above, the wakes are expected to have the effect of adding 
additional contributions of similar magnitude to each quadrant, therefore not contributing to an 
increase in the average Reynolds shear stress. The plots in figure [19] also show the fractional 
contributions from each quadrant when applying a threshold to the magnitude of the correlation 
u'fv'j. At both wall-distances the contributions from the second and third quadrants become 
completely dominant for thresholds above five times the product of the r.m.s. intensities of u'j and 
Wj. This reflects the large value of the negative skewness of the streamwise velocity pdf, i.e. the 
importance of the wake structures. 

The joint pdfs for the streamwise and wall-normal components of the particle velocity fluctu- 
ations at two corresponding wall- distances {y^ — 17 and 73) are shown in figure ^Oj The overall 
shape of the contourlines is found to be similar to the fluid counterparts in figure |18] yet events 
with highly negative values of the streamwise component are absent here. The fractional con- 
tributions from each quadrant to the correlation u'^v'^ (figure omitted) agree with the respective 
contributions of the fiuid velocity to the Reynolds shear stress. However, a faster decay of all 
quadrants' contribution with increasing threshold is observed for the particle velocity correlation 
u'^v'^, in accordance with the near-Gaussian shape of the pdf for the streamwise component (cf. 
figure 17 1). 



4 Lagrangian statistics 

In order to get further insight into the particle response to fluid motion, we consider statistics 
which describe the evolution of the particle velocity along individual particle trajectories. The 
Lagrangian auto-correlation of the particle velocity components as a function of the temporal 



separation r can be defined as follows (similar to Ahmed and Elghobashi 2001): 



where Uc,q(0 — ui%{t) — {uc){yj\xl!\{t) € Ij)Sai is the instantaneous velocity of the ith particle 
in the a direction, from which the mean streamwise particle velocity at the corresponding bin yj 
has been subtracted. By way of this definition, the correlation of the fluctuating particle motion 
can be studied while the contribution from the mean particle drift is eliminated. 



Our present results are shown in figure 21 It can be observed that all Lagrangian two-time 
correlation functions of case B exhibit a damped oscillation with a period of approximately 8h/ub 
superposed upon the typical decaying curves. For case A this feature is less pronounced, but it 
can still be observed with a period of approximately Ah/uf,, in particular in the wall-normal and 
the spanwise components. This phenomenon is believed to be related to the finite streamwise 
length of the computational domain. Considering that the difference between the mean velocities 
of the two phases is approximately equal to Uf, across most of the channel width (cf. figure [Tof )) , 
particles will encounter the same material fluid elements upon average after a relative turnover 
period of Trei — Lx/uh- Substituting the domain lengths of cases A and B (Lj, = 4/i and 8/i) 
yields a value for Trei which indeed matches the observed period separating the successive peaks 
of the Lagrangian two-time velocity correlation functions. However, these relative maxima of the 
correlation functions can only be caused by a significant interaction of particles with structures 
which are temporally coherent over time scales comparable to the turnover period T^ei- |Fede,| 



Simonin, and Villedieu (2007) have considered the problem of statistical bias caused by self- 
correlation of the fluid velocity along the path of heavy point-particles in homogeneous-isotropic 
turbulence (by means of one-way coupled simulations). They showed that choosing the parameters 
such that Trei > 4t£;, where te is the Eulerian integral time scale of the fluid velocity field, is 
sufficient in order to avoid bias of the Lagrangian correlation of the fiuid velocity 'seen' by the 
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particles. With the purpose of applying this criterion to the present case, the Eulerian integral time 
scale can be estimated from the Eulerian integral length scale and invoking Taylor's hypothesis. 
The streamwise integral length scale of the streamwise velocity component is defined as: 



L 



f-u 



(8) 



where i?„„ is the two-point auto-correlation function of fluid velocity. In single-phase channel 
flow, the data of Kim et al. ( 19871 yields Lf^u = 0.87h at the centerline {Lf^u = l-2/i at = 21). 
Supposing a convection velocity equal to the local mean velocity (Kim and Hussain 1993[ ), this 
leads to an estimate of te = 0-74h/ub at the center and 1.5h/ub in the buffer layer. Therefore, 
based upon the single-phase data the streamwise domain size in case B verifies the condition of 



Fede et al. (2007), while the domain in case A seems to be slightly too short in order to guarantee 
unbiased Lagrangian statistics. However, this a priori estimation does not seem to be able to 
explain that the bias effect is stronger in the case with the longer domain, where several repeated 
peaks up to the fifth one can be observed. Since a repeated positive correlation over four relative 
turnover periods requires that particles interact with fluid features of extreme longevity, the data 
shown in figure [21] suggests that there exist some large-scale coherent structures in the particulate 
flow which are not found in single-phase turbulence. In § [6] below we will indeed identify flow 
structures correlated over the whole box-length. 

Let us return to the discussion of the Lagrangian particle velocity correlations themselves. 
In the absence of reference data for channel flow, we will refer to the results of [Ahmed and| 
Elghobashi (2001) obtained by DNS of homogeneous shear flow with point-particles. Their most 
relevant conclusions with respect to the present work are: (a) particle inertia tends to increase 
the Lagrangian velocity correlations; (b) mean relative velocity (caused by gravity) leads to a 
decrease of the correlations due to the so-called 'crossing trajectories effect' ( Yudine[ 19591; (c) 
the anisotropic structure of the shear flow is reflected in the particle velocity correlations. 

Concerning the results shown in figure |21[ we remark that for all velocity components the 
decay is faster in case B than in the high-inertia case A, in accordance with the results of |Ahmed| 
and Elghobashi (2001 1. A rough match of the curves for both cases can be obtained by rescaling 



the temporal separations in case A with a factor of 1.5. 



We further observe from figure 21 that the streamwise (vertical) particle velocity component 



exhibits a much longer decorrelation time than the two cross-stream (horizontal) components, 
with a first zero-crossing occurring after as much as 85 (45) bulk time units in case A (B). The 
anisotropy is also evidenced by the integral time scales computed from the Lagrangian correlation 
functions, viz. 



' Lp,a 



RLp.air) dr. 



(9) 



which are shown in table [H It is well-known that the time scales of fluid structures in channel 
flow are strongly anisotropic. Choi, Yeo, and Lee (20041 have computed Lagrangian velocity 



correlations of fluid elements in the absence of particles. Their data at the center of the channel 
suggests a ratio of decorrelation time scales of roughly 1 : 0.4 : 0.5 for the three fluid velocity 
components w/ : Vf '■ Wf. By comparison, the corresponding correlations for the particle velocity in 
the present cases exhibit a much higher anisotropy of the streamwise component. This difference is 
believed to be related to the modified flow structure in the particulate case, where we have already 
observed a strongly enhanced intensity of the streamwise velocity fluctuations (cf. figure [TTj ). 
Unfortunately, no Lagrangian correlation data for fluid elements is available for the present cases. 
However, we will provide evidence of increased Eulerian two-point fluid velocity correlations in §[6] 
below. 

Concerning the crossing trajectories effect, the usual measure for gauging its importance is the 
ratio between the relative mean velocit y and the r.m.s. fluid velocity fluct uation s in the mean slip 
direction, 



\{uc) - {uf)\/^ (Csanady 



that the ratio for both cases ranges between 4 and 



19631. From figures 
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b) and 
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a) we deduce 



5 across most oTihe channel width. The 



14 



relative motion can therefore be judged as sufficiently strong for the particles to experience the 
consequences of their continuous change of the surrounding fluid environment, which is expected 
to contribute to a reduction of the temporal auto-correlation of particle velocities. However, in 
order to isolate the impact of the crossing trajectories effect upon the correlation functions, results 
from cases with different drift velocities and otherwise identical parameters would be necessary. 

The definition of the Lagrangian velocity correlation discussed so far does not distinguish the 
particles according to their location along the inhomogeneous direction. This distinction can be 
accomplished by modifying relation ([t]) such that the sums are performed only over those particles 
which are located within a wall-normal interval Ij at the initial time t. In doing so, one obtains 
correlation functions which turn out to have a similar shape to the globally averaged functions 



shown in figure^Tl but where the decay depends upon the particle's initial wall-distance. Figure 22 



shows the dependency between the integral time scale and the initial particle coordinate y in case 
B. The Lagrangian integral time scales associated with the streamwise and wall- normal particle 
velocity can be seen to increase progressively with the wall-distance; the spanwise component 
exhibits a similar behavior up to y ~ 0.8, then it decreases again towards the center of the channel. 
The trend of faster decorrelation near the wall is consistent with the fact that the length scales 
of the dominant turbulent structures in channel ffow increase with the wall-distance. IChoi et al.l 
( 2004 ) have shown for single-phase flow that this spatial inhomogeneity leads to a faster decay of 
the velocity correlations of fluid elements when they are released from smaller wall-distances. It 
is believed that a similar reasoning applies to our solid particles, with smaller fluid scales near the 
wall leading to a faster decorrelation of the particle velocity in that region. 



5 Spatial structure of the dispersed phase 

In this section we are concerned with the possible formation of instantaneous particle agglomera- 
tions in the flow. As discussed in the introduction, at least two mechanisms for clustering appear 
to be available: interaction with turbulent flow structures (i.e. 'preferential concentration') and 
the particle- wake effect. In order to determine whether agglomerations play a significant role, an 
analysis of inter-particle distances and local number densities is required. For this purpose we 
have applied four different diagnostic methods to our data: the average distance to the nearest 
neighbor, cluster detection, the probability of local particle concentration and the particle-pair 
correlation function. 

Let us define the average distance to the nearest particle as follows: 
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(10) 



where di 



is the distance between the centers of particles i and j. The time 



horn 



normalized by its value for a homogeneous distribution with the 
(||r2||/iVp)^/^, has been plotted in figure 
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Here the ratio 



evolution of the quantity d„ 
same solid volume fraction, d^. 

dmin/d!^^ has a lower bound of 0.2 (corresponding to each particle being in contact with at 
least one other particle) and an upper bound of unity (homogeneous distribution). Actual values 
recorded during our simulations vary mildly between 0.5 and 0.6, which is in fact very close to 
the value of a random distribution (numerically determined as 0.5621). It is interesting to note 
that direct simulations of pure sedimentation in triply periodic domains also yield values for the 
mean minimum inter-particle distance which (when normalized with dj^"™) fluctuate around 0.55 



(Kajishima 20041 



An alternative way of characterizing the spatial structure of the dispersed phase is by searching 
directly for the presence of particle clusters. Here we define a cluster as a set of particles of which 



each member is within a distance Ic of at least one other member (Wylie and Koch 2000). Since 



cluster detection is a computationally expensive task, it is necessary to resort to a fast numerical 



algorithm. We have implemented the technique suggested by Melheim ( 2005 ) which was applied 



to 400 (50) instantaneous particle distributions of case A (B), from which the probability of 
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the occurrence of clusters with different numbers of members n,. was determined for each case. 
Figure [24{^q, b) shows our results for two different cut-off lengths Ic = 2. 513 and AD. It can be 
observed that the pdf for the occurrence of clusters is a rapidly decaying function of the number of 
its members, and that the curve flattens as expected when increasing the cut-off length. Both DNS 
cases yield fairly similar results. The plots also include the corresponding pdf for a random particle 
distribution with the same solid volume fraction (numerically determined from 100 realizations), 
which is found to collapse with the DNS data. 

Since particles can interact over long distances via their wakes, the cluster definition was then 
modified in order to take the non-isotropic structure of the wakes into account. Instead of using a 
spherical cut-off radius, we define clusters through an ellipsoidal definition. More specifically, two 
particles i and j are considered as members of a cluster if their positions obey 



E 



''c,a 



< 1 



(11) 



for a given cut-off axis vector with components lc,a ■ Figure 24 [c) shows the pdf of the occurrence 
of clusters using a streamwise-elongated elliptical definition with a 3 : 1 : 1 axis ratio. It is found 
that the probability for large clusters indeed increases with respect to a spherical definition with 
the radius corresponding to the minor axis length. However, the random particle distribution 
again exhibits the same pdf as the DNS data to within statistical uncertainty. 

Several authors have used the statistics of local particle concentration as a measure of the 



tendency for preferential concentration ( 


Squires and Eaton 1991| Wang and Maxey 1993 


Fessler, 


Kulick, and Eaton 


1994 Aliseda, Cartellier, Hainaux, and Lasheras 2002 


). The concentration is 



determined by subdividing the flow domain into smaller boxes (with linear dimensions H^, Hy, Hz)-, 
counting the number of particles in each box, Up^box, and then dividing by the box volume, i.e. 
C = np^box/iHxHyHz). This process is repeated for each field, and finally statistics are computed 
over the homogeneous directions and the number of snapshots. In the present case, the width 
of all boxes in the wall-normal direction is taken as one particle diameter {Hy = £)), therefore 
allowing for a distinction of the results with respect to the wall-distance; the box-size in the two 
wall-parallel directions was varied between = Hz = 5D and AOD. Figure [25ja) shows the 
pdf of the local particle concentration in case B at a wall-distance of y+ = 28 for a cross-stream 
(horizontal) box-dimension of H^ = Hz ~ 20D. The concentration C is normalized by the average 
concentration at the same wall-distance, Cq. It can be seen that the DNS data agrees very well 



with a Poisson distribution which characterizes the probability of a random process ( Squires and 



Eaton 1991). The difference between the pdf from the DNS data and the Poisson distribution 



can be quantified by summing the square of the difference between both probabilities ( Wang and 
Maxey[|1993|), viz. 

2 



Do 



Pf 



(12) 



where P{n) is the probability of finding n particles in a given box. The quantity D2 is shown in 
figure 25 [b) as a function of the linear box-dimension and for different wall-distances. We observe 
that the difference always remains below 10^'^. Again, this definition can be extended to the 
non-isotropic case by defining elongated boxes (i.e. H^ > Hz). The results (which arc not shown 
here) lead to the same conclusion: the local concentration closely follows a Poisson distribution. 

As an additional measure of the spatial structure of the dispersed phase, we have determined 
the pairwise particle distribution function. This quantity corresponds to the probability of finding 
a second particle at a certain distance from a given particle. In order to investigate the two 
homogeneous directions separately, we have defined a directional (as opposed to radial) distribution 
function for streamwise and spanwise separations. In both directions, the curves are found to be 
essentially constants (the figures are omitted), indicating that all possible inter-particle distances 
are roughly of equal probability, irrespective of the wall-distance. 
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Finally we can conclude from the different diagnostic techniques employed here that the in- 
stantaneous spatial distribution of particles in the present cases is not significantly different from 
the distribution generated by a random process. This point will be further discussed in § [7] 



6 Coherent flow structures 

We have come across indications of a significantly modified turbulence structure with respect to 
single-phase channel flow at several points in the previous discussion. Here we attempt to identify 
the coherent flow structures which are responsible for the observed differences in the various 
statistical quantities. 

Figure [26] shows isosurfaces of the fluctuations of the streamwise velocity component for an 
instantaneous flow field in case B. The isovalues {u' — ±3.6wt-) are chosen such that in the cor- 
responding single-phase flow almost exclusively buffer-layer streaks are visualized. In the present 
case, however, the picture is fundamentally different. First, we recognize in figure |26[ b, d) the 
wakes trailing individual particles as streamwise elongated surfaces with negative fluctuation val- 
ues (due to fluid upflow and negative particle buoyancy). These wakes can be distinguished more 



clearly in figure 27 where a zoom into the graphs of figure |26[6, d) is provided. At the present par- 
ticle Reynolds number, particle wakes have a characteristic streamwise extension of less than one 
channel half width and a cross-stream extension comparable to the particle diameter. In addition. 



we observe in figure 26 very large coherent structures with streamwise dimensions of the order 
of the length of the current domain (8/i) and cross-stream scales comparable to the channel half- 
width. Very large structures with both signs of the isovalue for u' are found, the negative-valued 
structures being more difficult to discern due to the cluttering of the image by the large number of 
particle wakes. Visualizations of successive snapshots reveal that these fiow structures evolve on 
a relatively large time scale. An example of such a temporal sequence is shown in figure [28] where 
we can identify structures over an interval of approximately 45 bulk time units while they slowly 
revolve around each other in the cross-stream plane. This temporal coherence of the large-scale 
structures is at the origin of the slow convergence of the Eulerian statistics in the present case. The 
lack of symmetry of the streamwise momentum balance with respect to the midplane (observed 
3.1 figure [sji, b as well as figure 33 1 of appendix ]A]) is due to the marginal sampling of these 



large scales over the present observation interval. Furthermore, the presence of fiow structures 
with a lifetime of tens of bulk time units explains the repeated peaks of the Lagrangian particle 
velocity auto-correlation function discussed in § Jl] (cf. figure 21). Finally, it should be remarked 



that both particle wakes and the very large structures contribute to the strong increase in the 



streamwise r.m.s. fluid velocity, leading to the overall turbulence enhancement (cf. figure 111. 

In order to investigate whether the particles respond to the large-scale perturbations of the 
streamwise fluid velocity, we have visualized the instantaneous particle motion. The particles 
have been classified according to their streamwise velocity fluctuation: those with a high positive 
fluctuation value (w^ > 0.211^), highly negative fluctuation (v!^ < — 0.2uf,) and small fluctuation 
(— 0.2ub < < 0.2uf,). The instantaneous locations of the three classes of particles are shown in 
figures [29] and ]30] at the same time as the flow field in figure [26] The high-speed particles of both 
signs are clearly found in regions with a streamwise elongated shape, whereas the particles with 
streamwise velocity fiuctuation amplitudes below the threshold are more or less equally distributed 
everywhere else. Comparing the particle distribution with the locations of the fluid structures 



(figure 29 versus 26), a strong spatial correlation between high- intensity velocity fluctuations 
of the two phases is indeed observed. It has been confirmed by visualizing time-series that the 
spatial distribution of particles with intense streamwise velocity fluctuations exhibits the same high 
temporal coherence as the fluid counterpart. In particular, animations show that the large-scale 
organization of the flow starts to be visible approximately 12 bulk time units after introducing the 
particles into the flow. Subsequently, it remains a prominent feature over the entire observation 
interval. 

In order to further conflrm the statistical significance of the observed instantaneous fiow struc- 
tures we have computed Eulerian two-point correlations of fiuid velocity (taking into account only 
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points in the actual fluid domain) from 12 instantaneous flow fields. The data for case B is shown 
in figure 31 (for streamwise separations) and figure 32 (spanwise separations). For all components 
and in both spatial directions the decay at small separations is found to be enhanced in the present 
case as compared to the single-phase reference data. This result reflects the presence of smaller 
scales in the vicinity of the particles' surfaces as well as in their wakes. For larger separations 
the correlation functions of the particulate case and the single-phase flow cross over. The most 
pronounced features are the flnite correlation values of the streamwise velocity component for the 
largest separations. These correlation values are positive (negative) for large streamwise (spanwise) 
separations. This result is fully consistent with the presence of the observed large-scale structures. 
As a further consequence of the modified flow structure, all streamwise correlation amplitudes 
increase with the wall-distance contrary to what is found in single phase flow (figure 31). It is 
also remarkable here that for intermediate spanwise separations in the range of 50 < rj < 200 
the correlation functions for the streamwise and spanwise velocity components do not exhibit 
the characteristic negative minimum which is directly linked to the average spacing of low- and 
high-speed buffer-layer streaks (Kim et al. 1987). We still observe a weaker local minimum of 



the wall- normal component at the smallest wall-distance (figure 32 :i) . The evidence from the fluid 
velocity correlation functions shows that the structures which we have observed instantaneously 
have a large signiflcance for the flow statistics, masking in part even the signatures of the usual 
near-wall turbulence structures. 

The equivalent to the Eulerian two-point velocity correlation for the dispersed phase is the 
so-called particle pair velocity correlation (cf. Reeks, Fabbro, and Soldati 2006). It was conflrmed 
that the spatial decay of the particle pair velocity correlations (figures omitted) exhibits the same 
characteristic features as the counterparts for the fluid velocity. Thereby we have statistically 
substantiated the correspondence between particle motion and large-scale fluid structures observed 
in instantaneous flow flelds. 



7 Summary and discussion 

We have conducted a DNS study of turbulent particulate channel flow in a doubly-periodic domain, 
considering flnite-size rigid particles with numerically resolved phase interfaces. The simulations 
were performed in the dilute regime, allowing for an approximate treatment of direct particle 
encounters. In our simulations we have considered spherically-shaped heavy particles whose buoy- 
ancy was adjusted such that the terminal velocity matches the bulk flow velocity. The fluid was 
driven upwards by a mean pressure gradient, thereby fluidizing particles in the center of the chan- 
nel while particles near the walls moved downwards upon average. For the present bulk Reynolds 
number of 2700 and particle diameter of 1/20 of the channel half-width (corresponding to approx- 
imately 11 wall units), this leads to a terminal particle Reynolds number of 136. Two different 
density ratios were considered, differing by a factor of 4.5. The corresponding Stokes numbers of 
the two types of particles were C(10) in the near-wall region and 0(1) in the outer flow. The case 
with lower particle inertia was run in a computational domain with a length of 8 and 4 channel 
half-widths in the two wall-parallel directions, involving 4096 particles and integration times of 
C(IOO) bulk flow time units. 

We have presented statistical data for the Eulerian one-point moments of the velocities of both 
phases, velocity probability density functions, Lagrangian particle velocity correlations, various 
measures for the spatial structure of the dispersed phase, flow visualizations and Eulerian two- 
point correlations. The main results can be summarized as follows. 

1. The presence of particles strongly affects the carrier flow. On the one hand the mean flow 
profile tends towards a concave shape due to the interplay between the modified Reynolds 
shear stress profile and the mean force exerted by the particles upon the fiuid. The latter 
contribution is directly linked to the non-homogeneity of the mean solid volume fraction 
which has the effect of bringing more streamwise momentum towards the wall. On the other 
hand the intensity of turbulence is strongly enhanced with respect to the single-phase flow 
at the same bulk Reynolds number. This enhancement is mostly due to an increase in the 
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streamwise velocity fluctuations, whereby the anisotropy of the Reynolds stress tensor is 
considerably modified. 



2. At the chosen parameter values there is no significant clustering of particles. By 'clustering' 
we understand deviations from the local mean solid volume fraction beyond the extent that 
can be expected from a random process. 

3. The two dominating flow features are particle wakes and very large jet-like structures. The 
wakes constitute regions with low streamwise velocity (due to mean upflow and negative 
particle buoyancy), causing a marked negative skewness of the corresponding pdf. The 
large-scale structures appear as regions with positive or negative streamwise velocity fluctu- 
ations extending over eight channel half widths (the entire domain length) in the streamwise 
direction and approximately one half width in the cross-stream directions. As a consequence, 
the Eulerian two-point correlation of the streamwise fluid velocity exhibits finite values at 
large streamwise and spanwise separations. Furthermore, it was observed that the particle 
motion is strongly affected by these large-scale flow structures, whose signatures were found 
in Lagrangian particle velocity correlations and Eulerian particle pair velocity correlations. 

4. The differences between the two flow cases were relatively small. The fact that the high- 
inertia case was computed in a smaller computational domain might account partly for the 
observed differences. An exception are the Lagrangian particle velocity correlations which 
clearly exhibit longer decorrelation times in the high-inertia case. 

The fact that no particle agglomerations were found while the formation of new large-scale flow 
structures was detected merits further discussion. With respect to particle agglomerations it must 
be concluded that the two known mechanisms for their formation are not effective in the present 
context. For the preferential accumulation effect to generate signiflcant segregation of particles 
by centrifugal ejection from vortical structures, the particles need to be small compared to the 
dominant (i.e. high vorticity) flow structures. The particles studied here are probably too large in 
order to accumulate outside of the high- vorticity regions. This is particularly true in the near-wall 
region, where the present particles' size is comparable to the typical cross-stream dimensions of 
the streamwise vortices. 

Concerning the wake attraction effect, i.e. the formation of clusters due to the decreased drag 
experienced by trailing particles, it can be inferred from our results that the particle Reynolds 
number in the present study is probably too low for this mechanism to be of importance. In the 



pure sedimentation case (Kajishima 2004) particle clusters start to appear at a particle Reynolds 



number of 200; in the case with grid-generated turbulent background flow ( Nishino and Matsushita 



2004) trains of particles were found for particle Reynolds numbers of 140. It should be noted that 
in both of those studies the suspension was even more dilute than here. In turbulent channel flow 
wake-induced clusters have to our knowledge not been reported in the literature. Therefore, it 
would be interesting to perform future simulations with higher slip velocities in order to establish 
whether the wake-attraction mechanism is able to actuate in channel flow and to determine a 
critical Reynolds number. 

Even in the absence of particle agglomerations (i.e. signiflcant perturbations of the solid volume 
fraction) the flow is found to exhibit large-scale organization in the form of streak-like velocity 
perturbations. It should be emphasized that this feature is clearly particle- induced since no corre- 
sponding flow structures exist in single-phase channel flow. Since these large-scale flow structures 
were found to have a substantial effect on various statistical quantities (two-point velocity corre- 
lations, Lagrangian particle velocity correlations, r.m.s. fluid velocity fluctuations) the question 
about their origin is of great interest. It is plausible to believe that a mechanism related to hy- 
drodynamic instability is responsible for the formation of the large-scale flow structures. In this 
scenario it can be imagined that the addition of particles to the turbulent channel flow triggers 
an intrinsic instability eventually leading to the observed velocity perturbations. However, at the 
present time there is no conflrmation of the existence of such an instability mechanism. In this con- 



text it should be remarked that Matas, Morris, and Guazzelli (2003) found a decrease of the critical 
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Reynolds number when adding neutrally-buoyant particles above a certain size (D > dpipe/65) 
to laminar pipe flow. These authors speculated that the lowered transition threshold is due to 
the particle-induced velocity fluctuations triggering subcritical transition. The influence of the 
presence of dispersed particles upon the stability characteristics of pipe or channel flow deserves 
further investigation in the future. 

One further question which naturally arises in relation with the emergence of flow scales of 
the order of the streamwise period of the computational domain is whether the finite box size 
might play a decisive role. In particular, the typical streamwise size of these structures cannot be 
precisely determined if no decorrelation within the fundamental period takes place. It should be 
mentioned that these largest scales are not clearly observed in the case with the smaller domain 
(case A). On the other hand, computing the flow in much larger domains than the one used 
in case B is a challenging task, since those simulations are already on the limit of what can 
be afforded with the present algorithm on the largest available hardware. Nevertheless, we have 
recently initiated a companion simulation with twice the streamwise period of case B and otherwise 
identical conditions, corresponding to a Eulerian grid with more than 2.6 • 10® nodes. The results 
of this extended simulation — when available — will hopefully shed light on the scaling of the large 
structures. 
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A The effect of evaluating Eulerian statistics over the com- 
posite flow fleld 

In order to assess the effect of computing the Eulerian one-point statistics of the carrier phase 
over the composite flow fleld (including solid and fluid nodes), we have applied this procedure to 
12 instantaneous flow fields of case B which were taken over an interval of approximately 55 bulk 
time units. For the purpose of comparison, we have computed the corresponding quantities by 
averaging only over nodes which are located in the fiuid domain, using the same set of flow fields. 
This 'fiuid-only' average shall be denoted by (•)/. Both results are compared in figure 33 It can 
be seen that a significant difference is only obtained for the streamwise normal stress component, 
whereas the other stress components as well the mean velocity are not noticeably affected. The 
maximum normalized difference between the results of the two methods, as defined by 

^^^^^^^ maxy\{u)f-{u)\ ^ (13) 
maXy{u) f 

(and similarly for the other quantities) is given in table [6j It can be seen that the difference has a 
maximum of 6.7% in the case of the streamwise velocity fluctuation, while it is signiflcantly smaller 
for all other quantities. In particular, the difference for the mean velocity is below 0.7%. The 
larger difference for the quantities related to the streamwise velocity components is due to the fact 
that the virtual fluid velocity inside the heavy particles in upward flow represents principally a low- 
velocity perturbation with respect to the streamwise velocity of the surrounding fluid. Therefore, 
averaging over the composite flow field yields an underestimation of the mean velocity and an 
overestimation of the r.m.s. value of the streamwise velocity fluctuations. Figure [34] shows the 
difference between the two methods of computing these two moments of the streamwise velocity, 
plotted as a function of the wall-normal coordinate. It can be seen that the difference is directly 
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proportional to the local solid volume fraction (computed from the corresponding 12 particle fields) 
with a coefficient of 0(1). In particular, the least-squares fit of the data yields: 

{{u)f-{u))/ub « -1.4682(0,) (14a) 
{{u'u')f-{u'u'})/ul « 1.3634 ((/.,), (14b) 

with standard deviations of — 2.1966 • 10^'* and (t„„ — 2.7847 • 10~^, respectively. 

The result of the preceding paragraph is consistent with the following estimation. Let us 
assume that particles maintain a velocity equal to their terminal settling velocity with respect 
to the surrounding fluid velocity at all times and particle rotation can be neglected, i.e. the 
streamwise component of the virtual fluid velocity inside the solid regions can be expressed as 
m(x g iS) ss (^)/ ~ I'^c.ooh fluctuations of the solid volume fraction are also neglected. Then it 
can be shown that for the mean velocity: (u) — (u) f = — |itc.co|(0s). Likewise, for the streamwise 
component of the stress tensor: {u'u') — {u'u')f = —{u'u')f{(f)s) + \uc,ao\^{{4>s) — (<?^s)^); since we 
have imposed |itc,oo| = Ub (and the following two relations hold: (0s) <C 1 and {u'u') / <C Mb), the 
dominant term is |mc,ooP(0s). This suggests that for both quantities the difference scales with the 
local solid volume fraction, as was indeed observed in figure [34] 

From this comparison we can conclude that the statistics accumulated at run-time of the 
simulation, while not distinguishing between fiuid and solid locations, overestimate the r.m.s. 
value of the streamwise velocity fluctuations by up to 6.7%, while affecting the other fluctuating 
components much less and to an even lesser extent the mean velocity. The excellent fit of the 
present data with a linear function of the mean solid volume fraction enables us to correct the 
statistics of {u'u') which were accumulated at run-time. 



B The mean streamwise momentum balance 



Since the driving pressure gradient is not fixed in our simulations, we need to determine it first. 
For this purpose we consider the streamwise component of the momentum equation. 



dtu + (u • V)u + Vp 



(15) 



where u is the vector of fluid velocities, p the pressure normalized with the fluid density (p/), v 
the kinematic viscosity and f the volume force term which serves to impose the rigid body motion 
upon the fluid at the location of the particles. Integrating (15) over the whole spatial domain 



il and a sufficiently long temporal interval T to allow for a statistically steady state where the 
velocity field only depends upon the wall-normal coordinate, viz. 



dx 



2h + 2ul 



1 



ta+T 



to+T 



Np 

E 



fx dx dt , 



(16) 



where the last integral extends over the subvolume 5*^*^ occupied by the ith particle. The force 
integral can be eliminated by making use of the Newton equation for the linear particle motion 



(Uhlmann 2005a equ. 13), which states: 



Pp 
Pf 



1 u 



(m) 



f dx- 



Pp 
Pf 



1 Kg, 



(17) 



where Uc™^ is the vector of the mth particle's translational velocity. Substituting the streamwise 



component of (17) into (16) and making use of the fact that the fiow is statistically stationary 



(i.e. the average particle acceleration vanishes) yields the mean streamwise pressure gradient: 



dx 



' h 



Pp 
Pf 



gx- 



(18) 
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In order to derive the j/-dependent strearawise momentum balance, we now integrate equation ( 15 1 
over wall-parallel planes and time: 



where the Reynolds decomposition u = ((m),0,0) + u' has been substituted. Equation (19) could 
be used for evaluating the convergence of the statistics. However, the local average force density 
(fx) at the Eulerian grid nodes is not directly available in our data-sets. Alternatively, the quantity 



(fx) can be expressed by averaging the equation for the linear particle acceleration ( 17) over finite 
wall-normal bins, similar to the operation defined in Q. This procedure amounts to associating 
the sum of the forces applied to each particle with the particle's center location, i.e. the force 
density is computed as follows: 

(/.>(y.)« r r ir l E E / (20) 



The same operator applied to the equation of particle acceleration (17 1, assuming statistical sta- 
tionarity, yields: 

{fM-{j-^-^)9A<t>s){v,). (21) 



Figure 35 shows that the equality \2\\ is indeed verified across the entire channel by the data 



accumulated during our simulations. Finally, we substitute (21) into (19) to obtain the desired 
streamwise momentum balance equation ([5]) given in the main text (cf. § 3.1 1. 
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case 



Pf 

A 10 
B 2.21 



\g\h/ul St 



1.625 
12.108 



67 
15 



3.75 
0.83 



512 
4096 



n 

Ahx2hxh 
8hx2hx4:h 



Table 1: Physical parameters used in the simulation of particulate flow in a vertical plane channel. 
In all cases the particle diameter is chosen as £> = h/20, the global solid volume fraction is set to 
= 0.0042. 



case 

A 
B 



X Ny X 



1024 X 513 X 256 
2048 X 513 X 1024 



Atuii/Ax 

0.5761 
0.4148 



tobsUb/h 

700 
105 



Table 2: Numerical parameters employed in the simulations. Ni is the number of grid nodes in 
the ith direction, tabs is the observation interval after discarding the initial transient. The grid 



spacing in all cases is fixed at Ax 
per particle. 



/i/256, corresponding to Nl = 515 Lagrangian force points 



2/+ = 20 



present results, case B 
S{u'j:) Siv'j) S{w'f) 

-0.948 -0.135 0.019 
-1.289 0.029 0.031 



single-phase ( Moser et al. 1999 ) 

-0.323 -0.204 -0.006 
-0.629 0.030 -0.046 



Table 3: Skewness of the pdfs of fluid velocity components at two different wall distances, compared 



to single-phase reference data from Moser et al. (1999). 









5K) 


= 17 


0.089 


-1.221 


0.009 


= 73 


0.038 


0.003 


0.022 


= 220 


0.025 


0.014 


0.032 



Table 4: Skewness of the pdfs of particle velocity components in case B at different wall distances. 



case 

A 
B 



T, 



Lp,iUh/h 

13.58 
8.15 



Lp.2Ub/h 

0.86 
0.48 



T, 



Lp^sUb/h 

2.30 
0.79 



Table 5: Integral time scales (in bulk units) computed from the Lagrangian particle velocity 
auto-correlations, as defined in equation Q. 



£{{u)) £{^W)) SiVW)) SiVW^) £{{n'v')) 
0.0069 0.0672 0.0213 0.0070 0.0200 



Table 6: The maximum relative difference between the computation of the Eulerian statistics by 
averaging over the composite flow field and averaging over the fluid nodes only, estimated from 12 



instantaneous flow fields. The definition of the 'error' is given in ( 13 ) 
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Figure 1: The Cartesian grid in the vicinity of a sphere with diameter Dj i^x = 12.8, also showing 
the Lagrangian force points attached to its surface [Nl = 515). 




Figure 2: The geometrical configuration for vertical plane channel flow with suspended parti- 
cles. The fluid motion is induced by a negative mean pressure gradient in the direction of the x 
coordinate. The computational domain is periodic in x and z. 
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Figure 3: (a) The average vcrtieal partiele veloeity during start-up of the two-phase simulation. 
(6) The time evolution of the friction- velocity-based Reynolds number around the start-up time. 
Particles are added at t = to. case A; case B; single-phase flow. 
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Figure 4: The effect of discrete binning upon the data for the mean solid volume fraction (0s) in 

case B. (o) Profiles accumulated for different numbers of bins: Nbin = 40; Nbin = 80; 

Nbin = 160; - • Nbin = 320. (&) The data for Nbin = 320 in a close-up near the wall; the 
abscissa is normalized with the particle diameter. 
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Figure 5: The integral of the averaged streamwise momentum equation (|6]), evaluated from the 
statistical data accumulated over the observation interval in cases A (a and c) and B {b,d). (a), (6) 
show the full channel width; in (c),(d) the data is averaged over both channel half- widths, invoking 
odd symmetry. Statistics were accumulated over the composite flow field. The lines correspond 
to: Ttot + I4,; Ttou — o — y/h -I. 
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Figure 6: Wall-normal profiles of the raean fluid velocity, (a) in outer units, (b) in wall units: 

case A; case B; single phase flow ( Kim et aTj 1987). Statistics were accumulated 

over the composite flow field. In (b) — o — indicates the result of the single-phase computation 
with our code using the domain size of case A. 



32 




Figure 7: The wall-normal profile of the mean solid volume fraction for case A (- 
B ( ), using Nun = 80. 
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Figure 8: Wall-normal profiles of the derivative of the wall-normal particle velocity correlation; 
case A; case B. 
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Figure 9: (a) Wall-normal profiles of the Reynolds shear stress; line styles as in figure |6] In (b) The 
sum of the Reynolds shear stress and the particle contribution to equation ([6|, + J^, for 

case B ( ) is compared to the Reynolds shear stress of the single-phase flow ( ). Statistics 

were accumulated over the composite flow field. 
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Figure 10: Wall-normal profiles of (a) the mean particle velocity, and (6) the difference between the 

mean velocities of the two phases: case A; case B. Statistics for (u) were accumulated 

over the composite flow field. 
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Figure 11; Wall-normal profiles of the r.m.s. fluid velocity fluctuations; line styles as in figure |6j 
(a) Data from case A vs. single-phase flow; (6) case B vs. single-phase flow. The coordinate 
directions are indicated by the symbols: o a = 1; D a = 2; Aa = 3. The lines with the symbol 
"•" indicate the profiles of the streamwise component, corrected for the overestimation due to 
the accumulation of statistics over the composite flow field (cf. appendix |A]) ; errorbars show the 
standard deviation of the fit. 
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Figure 12: Wall-normal profiles of the turbulent kinetic energy k — X]?=i("i'"'i)/2; line styles as 
in figure |6] The lines with the symbol "•" indicate the profiles corrected for the overestimation 
of the streamwise stress component, due to the accumulation of statistics over the composite flow 
field (cf. appendix [a]); errorbars show the standard deviation of the fit. 
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Figure 13: Wall-normal profiles of the r.m.s. particle velocity fluctuations: (a) case A vs. case 
B with line styles as in figure [lOj (6) data of case B (solid lines) compared to the r.m.s. fluid 
velocity (dashed lines). The coordinate directions are indicated by the symbols: o a = l]'il\ a = 2\ 
A a = 3. In (6) the dashed line with the symbol "•" indicates the profile of the streamwise fluid 
velocity component, corrected for the overestimation due to the accumulation of statistics over 
the composite flow field (cf. appendix [A]); errorbars show the standard deviation of the fit. 
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Figure 15: Normalized probability density functions of fluid velocity components, computed from 
12 instantaneous flow fields of case B (taking into account only nodes inside the fluid domain) 
at — 20: (a) streamwise, (&) wall-normal, and (c) spanwise; present results (case B); 



single-phase data from Moser et al. (1999) at Rcr — 180; Gaussian reference curve. 
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Figure 17; Normalized probability density functions of particle velocity components in case B at 



different wall-distances: 
curve. 



17; 



■r 



73; 



y 



+ = 220; Gaussian reference 
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Figure 18: Isocontours of the joint probability density function of streamwise and wall-normal fluid 
velocity fluctuations in case B (evaluated from 12 instantaneous flow fields, taking into account 
only points in the fluid domain). The values for the contours vary between 0.0125 and 0.8 of the 
maximum probability, increasing by a factor of 2 between each level, (a) = 20; (6) y+ = 74. 
The dotted lines indicate \u'fv'f\ = —8{u'fv'f). 
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Figure 19: Fractional contribution from each quadrant to the Reynolds shear stress as a function 
of threshold: , first; , second; , third; , fourth, (o) y+ = 20; (6) = 74. 

The threshold is normalized by the r.m.s. intensities of each plane: ^{u'^u'j,) ■ ■^JWfV'f)- 
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Figure 20: Isocontours of the joint probability density function of streamwise and wall- normal 
particle velocity fluctuations in case B (evaluated from 250 instantaneous particle distributions). 
The values for the contours vary between 0.0125 and 0.8 of the maximum, increasing by a factor 
of 2 between each level, (a) wall-normal bin centered at = 17; (b) = 73. The dotted lines 
indicate \u'^v'J = —8{u'^vl,). 
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Figure 22: Integral time scales (in bulk units) computed from the Lagrangian particle velocity 

auto-correlations as a function of the particle's wall-distance at the initial time, case B. 

0.1 • TLp,iUb/h; TLp^2Ub/h; TLp^sUb/h. Please note the scaling factor applied to the 

data of the streamwise component. 
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Figure 23: Time evolution of the average distance to the nearest neighbor, normahzed by the 
value for a homogeneous distribution. Symbols as in figure fO ti marks the beginning of the 
observation interval. The dashed line corresponds to the value for a random distribution with 
the same solid volume fraction (computed for 524288 particles and correspondingly streamwise 
extended domain). 
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Figure 24: Probability density function of the occurrence of particle clusters with ric members 
in cases A (o) and B (•). The graphs correspond to different cut-off lengths: (a) spherical cut- 
off Ic — 2.5D; (6) spherical cut-off = 4Z); (c) elliptical cut-off with axes of length l^x = 7.5D, 
Icy = Icz = 2.5D. The probability was evaluated from 400 (50) instantaneous particle distributions 
in case A (B), spanning the respective observation intervals. The dashed lines correspond to the 
probability for a random distribution of particles with the same solid volume fraction, evaluated 
from 100 realizations. 
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Figure 25: Local concentration of particles in case B. (a) Probability density function of the 
number of particles per box of size 201? X D X 20D at a wall distance of y+ = 28: DNS results (•) 

versus Poisson distribution ( ); Cq is the average number of particles per box in the current 

y-slab. (6) The deviation between the pdf of the DNS data and the Poisson distribution according 



and z coordinate): 



to the definition (12), plotted as a function of the linear box dimension (in the directions of the x 
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28; 



y 



73: 



y 



220. 



50 



(a) (b) 




z/h z/h 



Figure 26: Instantaneous three-dimensional isosurfaces of streamwise velocity fluctuations u' — 
±S.6ur (equivalent to ±0.3uf,) of case B. The graphs (a, c) show surfaces corresponding to positive 
values, {b, d) show negative- valued surfaces. The view in (a, b) is into the wall, in (c, d) it is in 
the mean flow direction. Please note that only every eighth grid point in each direction was used. 
Here and in figures [28j [29j |30] the full computational domain is shown. 
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Figure 27: The same negative- valued isosurfaces as in figure [26[ b, d), but shown in the sub-domain 
[0,h] X [0,0. 5/i] X [0,h]. In addition, the particles are indicated by black spheres. 
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Figure 28: Isosurfaces of u' ~ ±3.6ut- in the cross-stream plane for successive snapshots with 
time increasing from top to bottom by increments of approximately 14.9ft./uf,. Positive-valued 
surfaces are shown in (a), negative ones in (b). The plots in the last row arc identical to those in 
figure [26]^ c, d). 
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Figure 29: Instantaneous positions of particles sorted according to their streamwise velocity fluc- 



tuation: (a) u'^ > +0.2ub; (b) u'^ < —0.2ub- The field corresponds to the one shown in figure 26 
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Figure 31: Two-point autocorrelations of fluid velocity fluctuations for streamwise separations, 
case B (evaluated from 12 instantaneous flow fields, taking into account only points in the fluid 
domain). The graphs show the following components of velocity: (a) streamwise, {&) wall-normal, 
(c) spanwise. The lines correspond to ?/+ ~ {21,87,225}, increasing along the arrows. The 



single-phase results of Kim et al. (1987) are indicated by the dotted lines 
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Figure 33: Eulorian statistics of case B, computed from 12 snapshots of the flow field by using 

two different methods: averaging over the composite field (including fluid and solid nodes); 

, averaging over the fluid nodes only, (a) The wall-normal profile of the mean velocity, (6) 

the r.m.s. velocity fluctuations, (c) the Reynolds shear stress. 
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Figure 34: The difference between the two procedures of accumulating the statistics shown in 
figure |33| (a) the mean velocity difference; (5) the streamwise normal stress component. The 
dashed lines indicate the proportionality with the profile of the solid volume fraction C • (<^s )(?/), 
where C ~ { — 1.4682, 1.3634} in (a) and (6), respectively, was obtained by a least-squares fit. 
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Figure 35; Wall-normal profiles of the mean streamwise force density (o) and the mean buoyancy 
term (x) for case B. 
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